rm(list = ls())

library(tidyverse)
library(lfe)
library(readxl)
library(broom)

## Helper functions

source("code/helper_functions.R")

## Load Data

btw <- readRDS("data/data_main.RDS") %>%
    dplyr::select(
        county_id_2018,
        year,
        matches("zweit|diff"),
    )


colnames(btw)

## Get treatment by remaining newspapers

treat <- read_rds("data/treat_to_n.rds") %>%
    dplyr::select(-state)

glimpse(treat)
treat$treat_categorical4 %>% table()

## Merge

btw <- left_join(btw, treat) %>%
    filter(!is.na(treat_categorical4))

## Reorder treatment

btw <- btw %>%
    mutate(treat_categorical4 = factor(treat_categorical4,
        levels = c(
            "NoChange",
            "DecAndInc",
            "Increase",
            "DecreaseTo1",
            "DecreaseTo2",
            "DecreaseToOther"
        )
    )) %>%
    filter(!year == 2013)


## Def outcomes

outcomes <- c(
    "small_party_zweit_diff",
    "polar_bund_alt_diff"
)


##

res <- lapply(outcomes, function(o) {
    ## Get DF

    df_temp <- btw

    ## Prep

    df_temp$outcome <- df_temp %>% pull(!!o)

    df_temp <- df_temp %>%
        filter(!is.na(outcome) & !is.na(treat_categorical4)) %>%
        dplyr::select(
            outcome, year,
            county_id_2018,
            treat_categorical4
        )


    ## Estimate

    m1 <- felm(outcome ~ treat_categorical4 | year | 0 | county_id_2018,
        data = df_temp
    ) %>%
        tidy_felm() %>%
        mutate(fe = "year")

    m1 %>%
        mutate(
            outcome = !!o
        )
}) %>%
    reduce(rbind) %>%
    mutate(period = "now") %>%
    filter(!str_detect(term, "Intercept|covar"))


## Remove some terms

res <- res %>%
    mutate(term = str_remove(term, "treat_categorical4")) %>%
    filter(!str_detect(term, "Intercept|covar"))

## Prep for plot


res <- res %>%
    mutate(outcome = str_remove(outcome, "_zweit")) %>%
    mutate(
        conf.low90 = estimate - qnorm(0.95) * std.error,
        conf.high90 = estimate + qnorm(0.95) * std.error
    )

## Rename

res <- res %>%
    mutate(outcome = dplyr::recode(outcome,
        `polar_bund_alt_diff` = "polar_bund_diff",
        `polar_land_alt_diff` = "polar_land_diff",
        `polar_unit_alt_diff` = "polar_unit_diff"
    ))

#### Main effects ####

res_plot <- res %>%
    filter(str_detect(term, "Decrease")) %>%
    mutate(outcome = dplyr::recode(outcome,
        `small_party_diff` = "Small party vote share",
        `polar_bund_diff` = "Polarization"
    )) %>%
    mutate(term = dplyr::recode(term,
        `DecreaseTo1` = "One outlet\nremaining",
        `DecreaseTo2` = "Two outlets\nremaining",
        `DecreaseToOther` = "Three or more\noutlets remaining"
    )) %>%
    mutate(term = factor(term, levels = unique(term)[c(1, 2, 3)]))


## PD

pd <- position_dodge(0.4)

## Plot

p1 <- ggplot(
    res_plot,
    aes(term, estimate)
) +
    geom_hline(yintercept = 0, linetype = "dotted") +
    geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90),
        position = pd, width = 0, linewidth = 1
    ) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
        position = pd, width = 0, linewidth = 0.5
    ) +
    geom_point(
        position = pd,
        shape = 21,
        fill = "white",
        size = 3
    ) +
    facet_wrap(~outcome, scales = "free_y", ncol = 2) +
    theme(legend.position = "bottom") +
    theme_bw() +
    xlab("") +
    ylab("Effect of newspaper exit") +
    theme(legend.position = "bottom")
p1
